The optimized Rayleigh-Ritz scheme for determining the 
quantum-mechanical spectrum 



Przemyslaw Koscik and Anna Okopiriska 
Institute of Physics, Swietokrzyska Academy 
Swi§tokrzyska 15, 25-406 Kielce, Poland 

The convergence of the Rayleigh-Ritz method with nonlinear parameters optimized 
bJQ' through minimization of the trace of the truncated matrix is demonstrated by a 

comparison with analytically known eigenstates of various quasi-solvable systems. 
■ We show that the basis of the harmonic oscillator eigenfunctions with optimized 

,— ,! frequency £1 enables determination of bound-state energies of one-dimensional os- 

Qh' dilators to an arbitrary accuracy, even in the case of highly anharmonic multi-well 

potentials. The same is true in the spherically symmetric case of V(r) = ^ + \r k , 

c3 ; 

if k > 0. For spiked oscillators with k < — 1, the basis of the pseudoharmonic oscilla- 
tor eigenfunctions with two parameters f2 and 7 is more suitable, and optimization 
J> ■ of the later appears crucial for a precise determination of the spectrum. 
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Since the early days of quantum mechanics, the techniques based on the variational 
j> ! principle have been successfully used for determination of ground-state characteristics of 
various physical systems. In its simplest form, the variational approach utilizes a single 
trial function that depends upon certain parameters, the values of which are fixed so as to 
minimize the expectation value of the Hamiltonian. This approach enables a very accurate 
determination of the ground-state energy, if the functional form of the trial function is 
appropriately chosen, but its extension to higher bound-states appears impractical because 
of difficulties in assuring orthogonality of the trial functions. A more suitable approach is 
offered by the linear Rayleigh-Ritz (RR) method with a trial function represented as a finite 
linear combination 

N-l 

<K^) = $> n Vv(T>), (l) 



n=0 
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where the functions ip n (~r*) are taken from some orthonormal basis in the function space. 
Treating the coefficients c n as variational parameters, one obtains a set of linear equations 

N-l 

Y,( H mn - s5 mn )c n = 0, m = 0, 1, N - 1 (2) 

n=0 

where H mn =< m\H\n > denotes the matrix element of the Hamiltonian between the states 
of the basis, \n >= ip n . The solution of the above equations yields iVth order approximations 
to wave functions </>^(~r^), and the corresponding energies e\ N \ for 2 = 0, N — 1 states. 
The accuracy of the RR method can be systematically improved by increasing the number 
N of basis functions, obtaining successive approximations to the larger and larger number of 
states. In this way a desired part of the spectrum may be determined with the approximate 
eigenvalues monotonically converging to the exact bound-state energies [l| at a rate strongly 
dependent on the choice of the basis. The method is also known under the name of exact 
diagonalization although the results become exact only in the limit as N tends to infinity. 
Long time ago, Hylleraas observed 2J that the effectiveness of the RR method may be further 
improved by introducing a scale parameter into the functions of the basis. Since then, various 
sets of functions depending on various nonlinear parameters have been used in determining 
spectra of atoms and molecules. The orthonormal sets made up of the eigenfunctions of a 
solvable Hamiltonian which depends on certain free parameters are especially convenient. 
In chemical applications the values of nonlinear parameters are fixed so as to a optimize 
the convergence of the RR estimate to the ground state energy; in the early works, the 
best values were found by trial and error [3( at present, they are fixed in computationally 
demanding procedures of iterative optimization jj, Q . 

Optimization of the RR scheme may performed on the basis of the principle of minimal 
sensitivity (PMS), which has been successfully used to improve perturbative calculation {(]- 
PMS requires the values of unphysical parameters introduced into calculation to be 
chosen so as to make approximations to physical quantities as less sensitive to the variation 
of these parameters as possible. This suggest to improve the RR method by fixing the val- 
ues of nonlinear parameters so as to minimize the Nth order approximation to the desired 
level energy E N . Such a strategy requires however diagonalization of the RR matrix to be 
performed in an algebraic way, which is feasible but only in low orders s|, or the appli- 
cation of extensive procedures of iterative minimization. Moreover, the scheme is not very 
economical, since the whole procedure must be repeated for each considered state. The op- 
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timized RR scheme, proposed by one of us [9j, adopts a different strategy that is also based 
on the PMS but insists on fixing the values of nonlinear parameters before diagonalization 
of the truncated matrix. Since the only physical quantity that can be determined before 
diagonalization of the Nth order RR matrix is its trace 

N-l 

Tr N H = < n\H\n >, (3) 

n=0 

which represents the sum of iV bound-state energies, we require the values of nonlinear 
parameters to be chosen so as to render Tr^H stationary The advantage of the scheme is 
that the iVth order approximations to many eigenstates are determined in one run and the 
obtained approximation to wave functions are orthogonal. It has been shown [9J that a good 
accuracy is acquired for the quartic anharmonic oscillator, using a basis of the harmonic 
oscillator (HO) eigenf unctions with optimized frequency. In this work we show that also in 
the case of other interaction potentials, the stationarity of the Tr^H condition optimizes the 
choice of nonlinear parameters. We note that implementation of the method with a modern 
software environment allows arbitrary precision calculations, which is especially important 
for determining the tiny energy splitting in the case of multi-well potentials with nearly 
degenerate levels. The computational cost of the optimized RR scheme is not high, as the 
50-digits accuracy of lowest bound-state energies is easily achieved with RR matrices of order 
N < 100, even in the most difficult cases of multi-well oscillators. Through the example of 
a sextic oscillator with analytically known eigenstates we show that the performance of the 
method for wave functions is also good, as various moments of the position operator appear 
well convergent. Later, we discuss the case of spherically symmetric potentials, where the 
two-parameter set of pseudo-harmonic oscillator (PHO) eigenfunctions constitutes a more 
appropriate basis for the optimized RR method. We gauge the convergence of the method 
for various anharmonic potentials Xr k , of both positive and negative power k, comparing the 
results with the exact solutions of different quasi-solvable examples (the sextic oscillator, the 
harmonium system, and the spiked oscillator). In all the cases studied, a good convergence 
is automatically ensured by using the trace condition for fixing the nonlinear parameters. 

The plan of our work is as follows. In section [Til we study the optimized RR method 
for one-dimensional anharmonic oscillators using the basis of the HO eigenfunctions. In 
section IHH the basis of the PHO eigenfunctions is introduced and the performance of the 
method for spherically symmetric potentials is discussed. Section [TV] is devoted to conclusion. 
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II. ONE-DIMENSIONAL CASE 



A. Harmonic oscillator basis 



The success of the RR method depends on the appropriate choice of the basis in the 
functional space. A basis constructed from the HO eigenfunctions 

= f^Lrl 2 H n (Qx)e-^ (4) 



with an arbitrary frequency Q playing a role of a nonlinear parameter h as p roved convenient 



for solving one-dimensional problems with purely discrete spectrum 
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-Il9j|. We use the 



above basis for determining spectrum of various anharmonic oscillators (AO). 



B. Quartic oscillator 

The most popular example is the quartic AO with the Hamiltonian operator given by 

where the units h = 1 and m = 1 are used. Bound-states of the system, tq hat exist if A > 0, 
correspond to those solutions of the Schrodinger equation 

H(f)(x) = E<f>(x) (6) 



which vanish at infinity. Mc Weeny and Coulson noted [10j that the HO wave functions 
(j3J) with any value Q > constitute an appropriate basis in the RR calculations for the 
quartic AO, as they fulfill the bound-state condition, but the convergence of the scheme 
depends strongly on the value of Q. They observed that the approximation series for the 
nth state energy converges quickly if the value of Q is fixed so as to minimize the matrix 
element < n\H\n >. However, such a prescription works well only in the case of a single- 
well AO (oo 2 > 0) but fails if the potential is double- well shaped (co 2 < 0) fll| . Moreover, 
diagonalization procedure has to be repeated with a different value of Q for each considered 
state. A more effective strategy is to take approximations to all the desired states from 
diagonalization of a single RR matrix with a compromise value of Q. Many ways of fixing 
the value of Q have been tested: so as to minimize the expectation value of the Hamiltonian 
in a chosen state of the HO basis (the first, namely |0 > in [ijj], the central in 13], the last, 



name 
in 



y \N > in [141 ]. that one for which the expectation of the Hamiltonian is the smallest 



Uj), so as to minimize a sum of expectation values in several HO states 15| . or so as 



to vanish the matrix element < 2|ff|0 > 16]. However, none of these prescriptions may be 
justified by the PMS, as none of the considered quantities does represent an approximation to 
a physical quantity. In this connection, one of us proposed jsj to fix the nonlinear parameters 
so as to make the trace of the truncated matrix stationary, since Tr^H represents the iVth 
order approximation to the sum of energies of the N lowest bound-states. For the HO 
basis PI, this amounts to fixing the frequency in iVth order calculation to the value ), 
which fulfils 



As shown in Ref. 



9| , the scheme automatically yields well-converging results for the spectrum 



of the quartic AO, in both the single- and double- well cases. We may add here that the values 
fiw , determined by the solution of ([7]), appear close to those for which the convergence of the 
above RR scheme is the quickest. This is demonstrated in Table HI where are compared 
with the values obtained by a numerical minimization of the error for the ground- 

state energy. A choice of nonlinear parameters has been considered on the same example 
of a quartic oscillator from a different perspective in Ref. [jij], where an analytic formula 
for optimum value has been derived from the asymptotic expansion. Comparing Table VIII 
of Ref. [l9| (where a corresponds to our Q 2 ) with our Table HJ we may see that the values 



derived from the asymptotic optimization scheme appear closer to than our values 
but an extension of this scheme to other systems and other basis sets is problematic, since 
asymptotic expansions which are uniformly valid in the nonlinear parameters are difficult to 
construct. Whereas, a fully algorithmic formulation of our optimization scheme allows an 
easy application for arbitrary systems. Generally, the accuracy of the optimized RR method 
diminishes with increasing number of the level, but the energies of N/2 lowest states can 
be trusted on, and using their values for calculating the free energy of the system, provides 



highly accurate results in a broad temperature range js]. Similarly, the well-determined 
part of the spectrum has been successfully utilized for approximate description of the time 
evolution in the quartic oscillator potential j^oj]. 

It is worthwhile noting that the optimized RR scheme may be implemented in a way 
allowing an arbitrary precision calculation by taking advantage of the present computer 
algebra abilities to deal with exact numbers. As an example we use the Mathematica 
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N 


y opz 


5E Q N (Qopt) 


y mm 


SEq ( ^771271 ) 


1 


1.29294233500847 


1.09716442402927 10 


-2 


1.29294233500847 


1.09716442402927 10~ 2 


5 


1.65920419620602 


8.63597923540916 10 


-7 


1.70670645005687 


1.33517756262170 10~ 7 


10 


1.86080663733626 


1.86447406929386 10" 


12 


1.95004181438340 


1.69012176051751 10~ 13 


15 


1.98940338014799 


6.71884044353837 10" 


19 


2.09841072018140 


2.7410652496 10~ 19 


20 


2.08593508969090 


4.36194667514212 10" 


-24 


2.20749144840388 


4.9101 10" 25 



TABLE I: The optimum value of the nonlinear parameter determined from the trace con- 

dition and the corresponding relative error of the ground state energy, 5E^ opt (Q opt ), compared 
with Q^A n , taken from Ref. 



for which the relative error of ground-state energy is minimal, 



package to calculate energy difference AE between the first excited and the ground state 
of the double-well oscillator with a Hamiltonian H = — + -^(x 2 — |) 2 for g = .001. 
Diagonalisation of the RR matrix of order N = 350 with the Mathematica precision constant 
set to 250, yields the value of the level splitting 



AE = 1.470464454175092501381989964494151981567800350052603 
52838605333378036605041575193505284182673433993282124674887 
3125608039469420665002677938176074660629119 10~ 68 , 



(8) 



which agrees with the 225 digits quoted with the result obtained on the basis of Zinn- Justin 
conjecture and confirmed by high-precision power series method 



Sextic oscillator bound-states 



The optimized RR scheme performs also well for oscillators with higher than quartic an- 
harmonicities, which present a more stringent test for approximation methods. We consider 
the sextic AO with a Hamiltonian given by 

72 , ,2 



c~ T I d 2 „ x , 



(9) 



that provides an interesting quasi-solvable example 22]. This is one of the rare cases when 



several bound-states of a system may be obtained exactly, namely p + 1 eigenstates are 
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known if the parameters of the sextic oscillator satisfy the condition 

u 2 = -(3 + 4p + 2i/)V2A, (10) 

where v = corresponds to the even-parity eigenstates and u — 1 to the odd-parity ones. 
The detailed explanation of how to obtain the exact eigenvalues and the explicit form of the 
corresponding wave-functions for a chosen value of p is given in the Appendix. 

Bound-state energies E(cu, A) depend on the parameters of the potential B 4|) . however 
in calculating the numerical results, we may set A = 1 without loss of generality. The 
eigenenergies at other values of A can be obtained as multiples of E(u\~^, 1), according 
to the scaling relation E(co,X) = A±E(u\~4,l), obtained by applying the transformation 
x i — y \~^x to the Schrodinger equation. As an example for testing the convergence of the 
optimized RR scheme we choose here the even parity case with p = 8, which according to 
( ITOl) corresponds to the quasi-solvable Hamiltonian 

h\ -- ld2 35v/ V + * 6 fin 

with a family of nine exactly known states (n = 0,2, ...,16). The iVth approximation to 
the AO bound-states, \n 3>, is obtained by numerical diagonalization of the Hamiltonian 
matrix in the basis of the N lowest even-parity wave functions of the HO (Tj0) with a fre- 
quency determined from ([7]) . In Table [TT1 and [TTT] the approximate values of energy 
and various moments of the position operator Xn^ N \^opt) =< C n\x k \n ^> are given for the 
ground-state (n = 0) and for the highest of the analytically known excited states (n = 16), 
respectively. A quick convergence of all the moments in Tables Hi] and III II indicates that 
wave functions are accurately determined. We have checked indeed that the agreement be- 
tween the approximate and exact wave functions to graphical accuracy is obtained already 
at N = 15. An exponential convergence for bound-state energies is evidenced in FigJU where 
the relative error 8 En (^opt) is plotted as function of N. In the case of the Hamiltonian 
H\ v= i = — — ^y^-x 2 + x 6 , when the family of nine odd-parity states is exactly known, 
the convergence properties of the optimized RR method are similarly good, but they slowly 
worsen with increasing value of u, when the double-well potential becomes deeper and more 
bound-states are known exactly. 

Also in case of higher anharmonicities, the RR method optimized by the trace condition 
enables highly accurate determination of the spectrum in relatively low order calculation. 
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For instance, diagonalisation of matrices of order N < 100 is sufficient for reobtaining all 



the bound states of the most difficult multi-well AO examples collected in [23| with the 
precision gradually diminishing from 50-digits for the ground state to 40-digits for the tenth 
excited state. The accuracy of various moments of the position operator is also good, as the 
hypervirial relations are fulfilled to a very high precision. Instead of presenting extensive 
tables of results, we refer the reader to our web page with the Mathematica program 24 ] 



which may be easily adapted for solving the desired AO problem to the desired precision. 



N 










20 
25 
30 
35 


-40.52625280948697 
-40.52625282342447 
-40.52625282343566 
-40.52625282343567 


2.72643166988877 
2.72643167388909 
2.72643167389268 
2.72643167389269 


23.60630747666730 
23.60630748254632 
23.60630748253900 
23.60630748253899 


242.23182746777999 
242.23182242437541 
242.23182241788116 
242.23182241787521 



TABLE II: The approximate energies EQ N \Q op t) and moments of the position operator x^ N \Q op t) 
for k = 2, 6, 10, obtained by means of the optimized RR method for the ground-state (n = 0) of 
the sextic oscillator (jlip at specific values of the dimension N. The underlined digits agree with 
the exact results. 



N 


EyPi^opt) 


X 1G \^opt) 






20 
25 
30 
35 
40 


40.52790329985854 
40.52625631425536 
40.52625282743947 
40.52625282343881 
40.52625282343567 


2.03702842810644 
2.03737443716647 
2.03737569324205 
2.03737569518445 
2.03737569518631 


35.34043326592744 
35.34277775827479 
35.34280112924913 
35.34280117889126 
35.34280117894960 


856.72311216403564 
856.42183863495585 
856.42125273736842 
856.42125207650909 
856.42125207596973 



TABLE III: Same as Table Hi] but for state n = 16. 



III. SPHERICALLY SYMMETRIC CASE 

The optimized RR method can be easily extended to the case of a central potential V(r), 
where r = |T^|, when the Schrodinger equation reduces to the one- variable problem at fixed 
angular momentum /. In the three-dimensional space, the problem is represented by the 
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FIG. 1: Semilogarithmic plot of the relative energy error, 5E^ (^opt)j f° r n=0 and n=16 state of 
the sextic oscillator (jlip in function of the dimension N of the optimized RR matrix. Here and in 
the following figures lines are drawn to guide the eye. 

equation 

\-i§T,r + + V(r)]R(r) = ER(r), (12) 

which upon introducing u(r) = rR(r) is transformed to the form 

H m u(r) = Eu{r) (13) 

with the radial Hamiltonian 



and boundary conditions u(r) — > 0, as r — > oo and u(0) = 0. 



A. Pseudoharmonic oscillator basis 



q The completeness of the set has been shown by Hall at al. |26j]. For A = 1(1 + 1), 
i.e. 7 = I + |, the solution (??) reduces to that of a spherically symmetric HO with a 
Hamiltonian 

^°--2^ + ^" + ^^- (15) 
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Optimization of nonlinear parameters of the PHO basis through the trace condition ([3]) 
amounts to choosing the values of Q and 7 in iVth order RR calculation so as to satisfy 

-^-Tr N H {l) = 0, and —TrjyH® = 0. (16) 
ait cry 

B. Radial oscillators 

Among interesting examples that can be easily treated by the optimized RR method are 
the radial oscillators described by the Hamiltonian of the form 

f>m 1 d2 1(1 + 1) u 2 r 2 , n . 

H w = — r + -i — H V Xr k (17) 

2rdr 2 2r 2 2 v ; 

with powers of anharmonicity k being both positive and negative. The above hermitian 
symmetric Hamiltonian is semi-bounded and consequently has self-adjoint extension in the 



Hilbert space L 2 (0, 00) which admits a spectral decomposition |27|]. Its spectrum is purely 

discrete. The Schrodinger equation for the above Hamiltonian may be rescaled in two 

1 2 
ways: the transformation r \— y r\( k+2 '> leads to the relation E(u, A) = A( fe + 2 ) E(z, 1) and the 

transformation r !->■ ru~^ yields E(u, A) = uiE(l, z ( 2 + ' ) with the dimensionless parameter 

2 

z = u\ ( fc + 2 ' . This shows that power k = —2 sets the border between very different behavior 
of the oscillators. Denoting by k a positive value | ju^g. |, in the case of k > —2 we have 
z = u\~ K and E(u, A) = \ K E(z, 1) and E(cj, A) = uE(l, z~ K ), which shows that for A — > 
or 00 0, the weak coupling or strong coupling perturbation theory may be respectively 
applied. The cases of k < —2 (spiked oscillators) are very different. For spiked oscillators, 
z = u\ K , E(u,X) = X K E(z,l) and E(u,X) = uE(l,z K ), which indicates that neither of 
the two terms of the interaction potential may be taken as dominant and the conventional 
perturbation expansions cannot be applied. The Klauder phenomenon [28] does not occur 
for radial problems with the boundary conditions of the Dirichlet type, i.e. u = at the 
singular point r = 0, and for A — > the energy eigenvalues of a spiked oscillator converge 
to those of the spherically symmetric HO of frequency 00 f lT5|) . However, the perturbation 
series contains the terms logarithmic in A and is ordered in fractional powers of A, thus 



unconventional methods have to be invoked for its derivation 



28 



29|. The perturbation 

theory with respect to Hho is supersingular because the matrix elements (m\r k \n) in the 
HO basis are infinite. Similar difficulties would be encountered for spiked oscillators in 



conventional RR calculations with the HO basis. As pointed out by Hall et si. 



30|, the 
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problem may be avoided by using the PHO basis with 7 > —k/2, since in this case the 
matrix elements appear finite and the RR calculation are well defined. For the anharmonic 
oscillators ( JT71) the numerical calculations may be simplified by removing the dependence 
on Q from the matrix elements by rescaling r h-> -^=, which yields 

= (2n + l)n5 mn + (m|r» - 

§(7 2 - 27 + I - i(« + l))H^|n) + ^Mr*|n>, (18) 
where the ket \n) corresponds to the radial function (??) with Q = 1, namely 

<(r) = {-l^J^^^^r-r-h-^F^ry). (19) 

In this work we consider three different examples: the spherically symmetric sextic os- 
cillator (k = 6), the harmonium potential (k = —1) and the "antisextic" spiked oscillator 
(k = —6), which all enjoy the nice feature of quasi-exact solvability We use the explicitly 
known solutions for testing the performance of the optimized RR method for both regular 
and singular spherically symmetric oscillators. In the case of positive power anharmonicities 
(k > 0), the optimum values of 7 turn out to be around I + 3/2, we may put therefore 
7 = 1 + 3/2 and use the basis of the radial HO eigenfunctions with frequency Q being the 
only parameter to be optimized. On the other hand, in the case of singular anharmonic 
oscillators (k < —1), the role of the parameter 7 is crucial, and we show that the values of 
'j opt determined from the trace optimization are always greater than —k/2, which ensures a 
successful calculation. 



1. Radial sextic oscillator 



First, we consider a sextic AO with the radial Hamiltonian of the form 

1 d 2 1(1 + 1) uj 2 o x fi 

As shown in the Appendix, the wave functions of the p + 1 lowest states are known in a 
closed-form if the parameters of the above oscillator are related as 



u 2 = -(5 + 4p + 2/)V2A, (21) 
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N 


P (JV), n \ 


2(JV) /n v 


6(iV)/ Q s 
'0 V s 'opt J 


10(JV)/ n \ 
'0 \} L opt) 


20 
25 
30 
35 


-48.11353417791531 
-48.11353418904174 
-48.11353418905196 


2.90220193255671 
2.90220193690299 
2.90220193690737 


27.98886643153595 
27.98886651687285 
27.98886651695577 


314.77282893841708 
314.77282672623041 
314.77282672244758 


-48.11353418905196 


2.90220193690738 


27.98886651695584 


314.77282672244337 



TABLE IV: The approximate energy Eq (Sl pt) and the moments of the radius operator 
rQ^ N \Q pt) (k = 2,6, 10) determined by the optimized RR method for the ground-state (n = 0) of 
the radial sextic oscillator (|22p for 1 = 1. 



N 






r 8* i^opt) 


^8 ^ ^ (ft opt) 


20 
25 
30 
35 
40 


48.11737731862358 
48.11354441726266 
48.11353420292910 
48.11353418906437 
48.11353418905197 


2.17521762647996 
2.17592836963907 
2.17593183660339 
2.17593184311651 
2.17593184312375 


42.02668954998873 
42.03132546920644 
42.03139851074727 
42.03139869962153 
42.03139869987743 


1136.40959750355092 
1135.67058214554528 
1135.66913552605609 
1135.66913419959634 
1135.66913419894809 



TABLE V: Same as Table IIV( but for the state n = 8. 



which is similar to that for the odd-partity one-dimensional case ( ITUj) but includes in addition 
the orbital number /. Besides a normalization factor the s-wave sector in a central 
potential is equivalent to the odd sector in a one-dimensional potential, therefore we need 
only consider the case I > 0. Setting A = 1, we discuss the case of p — 8, when 

m= ljjl + K^l) + (37 + 20V2 
2 dr l 2r A 2 

comparing the nine analytically determined eigenstates with the approximations obtained 
from diagonalization of the RR matrix in the basis of the radial HO eigenfunction with 
optimized frequency Q. The Tables IIVI and |V] contain our results for orbital number / = 1 
for the lowest (n = 0) and highest (n = 8) of exactly known states. It can be observed how 
the approximations to bound-state energies, El^\fl opt ) , and the moments r^^iflopt) =<< 
n\r h \n » for k = 2,6,10, approach the exact values with increasing N. In Figj2] the 
dependence of the relative energy error on A" is shown for both states. 
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FIG. 2: Semilogarithmic plot of the relative energy error, 5E^\Q opt ) for the two states (n=0,8) 
of the radial sextic oscillator (|22p as function of the dimension N of the optimized RR matrix. 

2. Harmonium 

The next example is related to the harmonium problem. Harmonium is a system of two 
particles confined in a harmonic potential and interacting via a Coulomb force, which enjoys 
the pleasant feature that the center of mass and the relative motion can be separated. The 
center of motion is subjected to a solvable harmonic oscillator equation, and the spherically 
symmetric relative motion equation corresponds to the k = — 1 power AO with a radial 
Hamiltonian of the form 

dr 2 r r 
The observation that the above Hamiltonian possess a closed-form solution 



31 



32| was 



of great importance and thus provides further rationale for investigation of approximation 
methods in the many-body theory. The exact solution exists if for some integer p, the 
frequency u fulfils the condition 

E = (3 + 21 + 2p)u, (24) 

which, in contrast to (I2TI) . does not depend on the coupling A but on the energy of the 
investigated state. In difference with the case of the sextic oscillator, for harmonium there 
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is but a single bound-state, not necessarily the ground-state, that is known exactly. This 
happens if the parameters u and A are related by a particular relation that depends on 
an integer p, as given in the Appendix (1A5jl . By utilizing the scaling property E(u,X) = 
X 2 E(u>\~ 2 , 1) and setting A = 1 we have that Eq. ( 1A5I) determines the values u p at which 
the bound-state is analytically known. Using the exact bound-state solutions for testing 
the convergence of the optimized RR method in the PHO basis, we compare three ways of 
fixing the nonlinear parameters: first, a naive RR scheme (Q = u p and 7 = 3/2 + 1), second, 
optimization of the parameter Q (Q = Q opt and 7 = 3/2 + /), and third, optimization of both 
parameters (Q = Q opt and 7 = ■y op t)- Figj3] presents the semilogarithmic plot of the relative 
energy error 5E^p as function of N for the largest possible value of confining frequency for 
which the ground state solution is known, namely u>i = 0.25 (I = 0,p = 1). A similar plot 
for a smaller value of frequency, co^ = 35 ~ 4 3 2 ^ ~ 0.00867 (I = 0,p = 4), is given in FigJU 
The results for higher values of angular momentum I show a similar tendency, therefore we 
do not show them. In the case of large confinement frequency the optimized scheme proves 
superior over the naive one, but in the strong correlation limit (a; << 1) it is the naive one 
that works slightly better. We may notice that the application of the trace condition in the 
RR scheme ensures an exponential convergence, albeit its rate is generally slower than in 
the case of positive power AOs. Since the values of nonlinear parameters determined from 
the trace condition appear close to 7 = 3/2 + I and Q = u p , the optimized RR scheme may 
be regarded as a justification for using the naive RR method for harmonium-like systems. 



Now, we come to the computationally more difficult case of a spiked oscillator with 
anharmonic potential of the power k < —2, which exhibits a highly singular behavior at 
the origin. As an example we consider the "antisextic" oscillator (k = —6) with the radial 
Hamiltonian given by 



A bound state of the system may be considered analytically determined if its energy is given 



3. Spiked oscillator 




(25) 



by 



E ip) = 2(p+l)w 



(26) 



15 




FIG. 3: Semilogarithmic plot of the relative energy error for the ground state of harmonium with 
confinement frequency u\ = 0.25 as function of the dimension N of the optimized RR matrix. 

and a specific relation between u and A (IA16I) is satisfied for an integer p. As opposed to 
the previously discussed sextic oscillator and harmonium cases, where the results have been 
rescaled in terms of A, in the case of spiked oscillators we exploit the rescaling in terms of u, 
namely E(u, A) = uE(l, u 2 X). This allows us to put ou — 1 in the numerical calculation in 
accordance with the usual practice in the literature on spiked oscillators. The specific values 
of A at which the bound states are analytically known may be obtained by substituting 
(1A17j) into (1A16I) and determining the p + 1 solutions of the polynomial equation ( 1A5[) . If 



A p (j) are numbered with decreasing value of A from i = to i = p, then the ith exactly 
known state represents the ith excitation in the respective potential. 

The PHO basis has been used in the RR calculations for spiked oscillators by Hall et 
al. {26], who derived useful formulas for the matrix elements of the operator r k in the 
PHO basis on condition of 7 > —k/2. In the case of the "antisextic" oscillator this 
amounts to 7 > 3, and the expressions for matrix elements (m\Wn) are singular at the 



values 7 = 0,1, 2, 3. Hall et al. [26] used the prescription of fixing nonlinear parameters 
(A = 7 2 — 27 + I and B = Q 2 ) so as to minimize the approximate energy of the con- 
sidered level -E^(f2, 7). In low orders, an algebraic diagonalization of the RR matrix has 
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FIG. 4: Like Figj3l but for confinement frequency u>4 ~ 0.00867. 

allowed them to obtain analytical approximations to ground-state energy of the " antisextic" 
oscillator, but for its precise determination in higher order calculation, a time consuming 
procedure of iterative optimization of nonlinear parameters must have been used by Naser, 
Hall and Katatbeh [5]. Our scheme requires much less computational cost, since the opti- 
mum values of nonlinear parameters are determined from the trace condition and further 
diagonalization of the RR matrix is performed only once in each order calculation. The 
convergence of our method is demonstrated on two quasi-solvable examples of the "antisex- 
tic" oscillator with angular momentum I = 0, where ground-states energies are analytically 
known. The semilogarithmic plot of relative errors is shown in Figj5] for p = 2 case ( 
A 2(0) = ^ [9887 + 32V333778 cos[| arctan( 18 S^ jP)]l « 369.26, E$ = 6), and in Figj6] 
for p = case (Ao(o) = y§8 ~ 0-07, Eqo = 2)- We observe that for both large and small 
values of A, the results of the two-parameter (D, op t, jopt) an d the one-parameter optimization 
(Q = u = l,7opt) are nearly the same. We have checked that for other spiked oscillators 
the case is similar and the calculations may be simplified, since only parameter 7 needs to 
be optimized. This is in difference with the results of the approach that utilizes iterative 
optimization of ground state energy, presented in Table II of Ref. jsj, where much quicker 
convergence has been obtained by optimizing both parameters. However, plotting the rel- 
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ative errors for exited states determined from diagonalization of the RR matrix with the 
values of parameters taken from Table II of Ref. 5| for the case of A = 0.1 in Fig. [7J we 
observe that the precision of energy determination in this approach rapidly decreases as the 
number of the level increases. On contrary, the results of our approach, obtained from the 
RR matrix of the same dimension N = 80, plotted on the same figure, indicate a uniformly 
good precision in a wide range of energy levels. We conclude that the values of nonlinear 
parameters obtained from minimization of the trace of the RR matrix are appropriate for 
a precise determination of the whole part of the spectrum, although this is not necessarily 
the best possible choice for a particular level. 
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FIG. 5: Semilogarithmic plot of the relative energy error for the ground-state of the radial spiked 
oscillator (|25p of the strength A2(o) ~ 369.26 as function of the dimension N of the optimized RR 
matrix. 

In Table IVTl and I VIII we show the numerical results of the RR method with the parameter 
7 optimized through minimization of the trace for the above- discussed quasi-solvable cases. 
One can observe how the values of ground state energy and various moments of the radial 
position operator converge to the exact values with increasing dimension of the RR matrix 
N. The optimum value of 7 depends strongly on A and grows with N. In all the cases we 
studied, the condition 7 op4 > 3 is fulfilled, although the smaller A is, the closer -y opt is to the 
value 3, where the matrix elements become singular, which explains the heavy worsening of 
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FIG. 6: Like FigEl but for the strength A 0(0 ) ~ 0.07. 

convergence with decreasing A. Nevertheless, we may observe that for A — >■ the bound- 
state energies approach smoothly those of the radial HO, which is due to using the PHO 
basis that ensures the satisfaction of the Dirichlet boundary condition in the RR calculation. 



N 


lopt 




Xq i (Tout) 


xl (N) {lopt) 


x 6 {N) (7ppt) 


20 
40 
80 
120 


14.48 
17.29 
20.97 
23.61 


6.00021390368223 
6.00000223509568 


2.84582121384655 
2.84562114647890 


8.28291890988045 
8.28186647605726 


737.684487549365 
737.249285331635 


6.00000006213529 


2.84561903526830 


8.28185556086959 


737.240961117416 


6.00000000480818 


2.84561898882020 


8.28185531621157 


737.240889729308 


exact 




6 


2.84561898466095 


8.28185529459909 


737.240860856683 



TABLE VI: The approximate energy Eq (j pt) and the moments of the radius operator r (7opt) 
(k = 1,2,6) determined in the optimized RR method for the ground-state of the radial spiked 
oscillator ([25]) of the strength A 2 (o) ~ 369.26. 
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FIG. 7: The semilogarithmic plot of the relative error of bound state energy of the spiked oscillator 

n 

of A = 0.1 as function of the level number n, determined by the method of Ref. [5J and by the 
optimized RR method of the present work. 



N 


lopt 




r { N) {7o P t) 




rl {N \lopt) 


20 


4.38 


2J) 1776622625948 


L44367633236718 


2.28981024534251 


30.6645428912254 


40 


5.02 


2.00734163831413 


1.43695960594035 


2J27129896763614 


29.8195939603567 


80 


5.88 


2.00272386999177 


1.43401033592011 


2^26321009528890 


29.6572447265661 


120 


6.50 


2.00140501517273 


1.43316756712602 


2^26090435812050 


29.7830590372660 


exact 




2 


1.43226578557733 


2.25844053161144 


29.4482015786915 



TABLE VII: Like Table EH but for the strength A 0(0) ra 0.07. 



4- Generalized oscillators 

We have tested the convergence properties of our approach further by considering the 
potential V(r) to be a linear combination of r s and r* with the power s being negative and 
t being positive. The optimized RR method does allow solution to eigenvalue problem for 
various combinations of potential parameters. In FigJH] we have plotted the error of ground 
state energy determination for the exemplary potential V(r) = r s + r* with various powers 
s = —1.5, —1.8, —1.95, —2.05, —2.2, —2.5 and t = 2, 4, 6 in function of the dimension of the 
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optimized RR matrix, N. The convergence becomes exponential at not too large values of N, 
with the rate depending on the detailed shape of the potential. In all the cases considered, 
using the PMS condition for the trace of the RR matrix for fixing the nonlinear parameters 
Q and 7 ensures an effective determination of the spectrum. 
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FIG. 8: The semilogarithmic plot of the relative error of the optimized RR method for bound state 
energies of the generalized oscillator V(r) = r s + r* for various powers s and t, as a function of the 
dimension N. 



IV. CONCLUSION 



We have discussed optimization of the RR scheme by introducing nonlinear parameters, 
where values are fixed by minimization of the trace of the truncated matrix. Using the basis 
of the HO eigenf unctions with optimized frequency Q, we obtain an efficient method for 
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determining spectrum of multi-well one- dimensional AOs to practically any precision. In 
the case of radial oscillators with \r k anharmonicity, the basis of the PHO eigenfunctions 
with two arbitrary parameters Q and 7 seems more suitable. For positive power oscillators 
(k > 0) the role of the parameter 7 turns out to be minor, and the scheme may be simplified 
by using the radial HO eigenfunctions (i.e. setting 7 = 3/2 + /). In the case of negative power 
oscillators k < — 1, the parameter Q plays a minor role and may be set equal to the frequency 
u of the harmonic term in the Hamiltonian, whereas optimization of the parameter 7 in each 
order calculation is crucial for a good convergence. In the limiting case of harmonium-like 
potential (k — — 1) the optimum values of both parameters turn out to be equal to the ones 
that are used in the naive RR method (Q = u and 7 = 3/2 + /). The optimized RR method 
performs well not only for energies but also for wave functions, yielding well-convergent 
approximations to various moments of the position operator. 

The RR method optimized by the trace condition appears effective for numerical calcu- 
lation and may be used to arbitrary accuracy within the modern software environment. It 
turns out that far greater improvement of accuracy is obtained from optimizing the nonlin- 
ear parameters in each order calculation by the trace condition rather than from increasing 
the dimension of the basis with the parameters remaining fixed. The computational cost of 
our scheme is much lower than in the case of iterative optimization of nonlinear parameters. 
Another advantage is that the whole set of energy levels may be determined at once and the 
approximate eigenvectors are mutually orthogonal. For the class of potentials with purely 
discrete spectrum considered in the present work, the RR method with the PHO basis is 
highly competitive with existing methods, the results of which are easily recovered in our 
approach. The convergence of our approach may be slower than achieved in the methods 
specialized for a particular problem, but the advantage is that the results are obtained auto- 
matically for a large class of systems by the algorithm described above without the necessity 
of specifying any starting value. 



Appendix A: Quasi-exact solutions 

In the appendix we consider the problem of quasi-solvability of anharmonic oscillators, 
thereby deriving the formulas for the analytically known solutions that are used for testing 
the convergence of the optimized RR method. Generally, the problem is quasi-solvable if 
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the eigenf unction can be represented as 

p 

V>(i/) = f(y)J2 a ™y n > 

n=0 

and the coefficients of the series satisfy the three-term recurrence relation 

A n a n+ i + B n a n + C n a n _i = 0, n — 0, 1, 2, ... 



(Al) 



(A2) 



where a_i = 0. The coefficients of the series a n as well those of the recurrence A n , B n , C n 
depend solely on the parameters of the potential and the bound-state energy E. In order 
for the series in (IA1I) to terminate after the p— th term, we must have 



and 



C»+i — 



a p+1 = 0. 



(A3) 



(A4) 



It is easy to convince oneself [33] that the second condition is equivalent to 
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(A5) 



\ C p B p J 

where the determinant is a polynomial of the degree p + 1 in the variable E. The two equa- 
tions ( 1A3I) and (1A5I) may be used to determine the specific relations between the parameters 
of the potential that must be satisfied for the exact solution to be of the form (1 All) . 



1. One-dimensional sextic oscillator 

For the one- dimensional sextic AO ( 1IIIB4I) the analytically known eigenfunction, in the 
even-parity [y = 0) and odd-parity (v = 1) case, may be represented as 

A 4 P 

^(x) = a n x 2n+ \ (A6) 

n=0 
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where the coefficients of the series satisfy the recurrence relation (1A2|) with 

A n = (2n + 2 + v)(2n + l + v), 
B n = 2E, 

C n = -co 2 - (4n - 1 + 2v)V2\. (A7) 

The series in (1A6I) terminates after the p— th term, if C p+ i = and a p+ i = 0. The first 
condition is fulfilled if the parameters of the sextic oscillator satisfy 



u 2 = -(3 + 4p + 2z/)V2A. (A8) 

For a chosen value of p, this corresponds to two quasi-solvable cases of a sextic oscillator: i) 
v = 0, where the p + 1 lowest even-parity states are known; and ii) v — 1, where the p + 1 
lowest odd-parity states are known. In both cases, we can obtain the p+ 1 exact bound-state 
energies E(u, A) by solving the polynomial equation f]A5[) . and thus the corresponding exact 
wave-functions by determining the non- vanishing coefficients a n from the recurrence relation 
2D. 

2. Radial sextic AO 

For the spherically symmetric sextic AO (|2"U|) the analytically known eigenfunction be- 
come 

A 4 P 

u{r) = r m e~^3 a n r 2n } (A9) 

n=0 

where the coefficients a n satisfy the relation f ]A2j) with recurrence coefficients 

A n = 2(n + 1)(3 + 2n + 2/), 

5 n = 2E, (A10) 
C n = -w 2 - (1 + 4n + 2/)\/2A. 

The closed-form solutions exist if the sextic oscillator parameters satisfy 

co 2 = -(5 + Ap + 2l)V2\ = 0, (All) 

and the p + 1 bound-state energies are determined by the condition (1A5|) with recurrence 



coefficients of the form flAlOl) 
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3. Harmonium 

The eigenfunction of the harmonium-like Hamiltonian (|23|) can be written as 

p 

u{r) = r m e"^ 2 G ™ r ™' ( A12 ) 
n=0 

where a n satisfy the relation (1A2j) with recurrence coefficients 

A n = (n + l)(n + 2l + 2), 

C n = £ - (1 + 2n + 2l)u. (A13) 

The closed- form solution is obtained if C p+ \ = 0, which means that 

E = (3 + 21 + 2p)uj, (A14) 

and the condition (IA5I) is satisfied with the recurrence coefficients given by flA13j) . Compared 
to the anharmonic oscillator case, the condition (IA14I) depends on energy; therefore for a 
particular value of u, denoted by u p , only the bound-state with energy E p = (1 + 2n + 2l)u p 
is known exactly. 



4. Spiked oscillator 

For the spiked AO with the radial Hamiltonian of the form (1251) . the analytically known 
eigenfunction assume the form 

2 P 

u(r) = r^e 2 \ a n r , (A15) 

n=0 

where the coefficients a n satisfy the relation (IA2[) with recurrence coefficients given by 



A n = -4V2A(n + l) 

B n = -- - 4n(l + n) + 1(1 + 1) + 2^0; 

C„ = -2(E-2nuj). (A16) 
Closed-form solutions are obtainable if 

E = 2(p+l)u, (A17) 
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and the condition f ]A5j) is satisfied with recurrence coefficients as stated in ( 1A16j) . For a 
fixed p, Eq. (IA5j) possess p + 1 solutions that determine the p + 1 cases of specific relations 
between A and u, at which the bound state of energy (1A17I) is analytically known. 
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